This repository contains tools for single cell simulations in barrel cortex.
installation instructions, overview of this package¶
Setup¶
If you are in the IBS-group, you can activate an environment, that already contains all of the packages mentioned below:
source /nas1/Data_arco/prgr/anaconda2_florida/bin/activate in_silico_frameworkTo get a fresh installation, follow the following steps:
- download and install Anaconda 2.7 https://www.continuum.io/Downloads
- install neuron such that you can import it as python module. Version 7.4 has been tested, newer versions are likely to work. Detailed installation instructions can be found at http://www.davison.webfactional.com/notes/installation-neuron-python/. If you use Ubuntu and you have trouble compling neuron, it has been reported that installing the following packages solves the problem:
sudo apt-get install bison flex g++ libxt-dev xorg-dev python-dev libncurses5-dev - Add the neuron folder to your PATH environment variable, such that you can run
nrnivmodlanywhere - install the following dependencies:
- sumatra, used for parameterfiles:
pip install sumatra - pandas 0.19.2, data analysis library:
conda install pandas==0.19.2 - dask 0.16.1, dynamic task scheduling and "big data" extension of pandas:
conda install dask==0.16.1#was 0.14.3 - distributed 1.20.1 allows non-bloccking computations and brings dask to a cluster:
conda install distributed=1.20.1#was 1.15.2 - seaborn: statistical data visualization:
conda install seaborn==0.8.0 - fasterners: robust file based locking:
pip install fasteners - jinja2: html template engine, required for embedded animations:
pip install jinja2
- sumatra, used for parameterfiles:
- clone or pull this repository:
git clone https://github.com/abast/in_silico_framework.git. - add your in_silico_framework folder to the PYTHONPATH variable
- unzip the following folder: in_silico_framework/getting_started/barrel_cortex.zip such that the following file exists:
in_silico_framework/getting_started/barrel_cortex/nrCells.csv - Open the file model_data_base/simrun2/seed_manager and adjust the variable
path. This will specify the location where used seeds are saved. Any location where you have write access is suitable. run the test suite:
python run_tests.py.Due to the statistical nature of the model, some tests might fail from time to time. These tests have the word statistical in their description. If such a test fails, run the testsuite again. If that test fails again, there most likely is an issue. Tests, that do not have a statistical flag in their description may never fail.
conda install -c anaconda lz4
conda install -c anaconda blosc
conda install -c conda-forge python-blosc
#conda install pytables
pip install mockconda install -c anaconda gcc_linux-64
conda install -c psi4 gcc-5
conda install -c conda-forge binutils
conda install -c conda-forge gitdask:0.14.3 distributed: 1.6.3 sklearn: 0.18.1 cloudpickle: 0.2.2 pandas: 0.19.2
dask:0.16.1 distributed: 1.20.1
Packages in this repository you will most likely directly interact with¶
singlecell_input_mapper: Module for generating anatomical models, i.e. determining number and location of synapses and number and location of presynaptic cells
single_cell_parser: High level interface to the NEURON simulator providing methods to perform single cell simulations with synaptic input. The anatomical constraints of the synaptic input are provided by the single_cell_input_mapper module.
simrun2: High level interface to the single_cell_parser module providing methods for common simulation tasks. It also provides methods for building reduced models mimicing the full compartmental model.
model_data_base: Flexible database whose API mimics a python dictionary. It provides efficient and scalable methods to store and access simulation results at a terrabyte scale. It also generates metadata, indicating when the data was put in the database and the exact version of the in_silico_framework that was used at this timepoint. Simulation results from the single_cell_parser module can be imported and converted to a high performance binary format. Afterwards the data is accessible using the pandas data analysis library and dask.
single_cell_analyzer: Library for analysis of single_cell_parser results. Use this module, if you specifically want to analyze a single simulation run. If you want to analyze the results of many simulation trails, the recomended way is to import the simulation results in a model_data_base and use pandas and dask for the analysis afterwards.
Other packages¶
mechanisms: NEURON mechanisms (Ion channels, synapses, ...) used by the neuron simulator. If you run import mechanisms, these mechanisms are beeing compiled. Make sure, that you have the nrnivmodl executable in your path, otherwise this will not work
test_simrun2, test_model_data_base, ...: tests for the respective modules. To run the testsuite, use python run_tests.py
How should I interact with these packages?¶
The recommended way is to use the Interface module which glues together all these packages mentioned above to one Application: It provides the API necessary to perform simulation tasks and provides additional methods that improve interactivity.
import Interface as I
Now, you can access the relevant packages, functions and classes directly:
I.scp # single_cell_parser package
I.sca # single_cell_analyzer package
I.ModelDataBase # main class of model_data_base
I.map_singlecell_inputs # compute anatomical model for a given cell morphology in barrel cortex
I.simrun_run_new_simulations # default function for running single cell simulations with well constrained synaptic input
I.mdb_init_simrun_general.init # default method to initialize a model data base with existing simulation results
I.mdb_init_simrun_general.optimize # converts the data to speed optimized compressed binary format
I.synapse_activation_binning_dask # parallelized binning of synapse activation data
I.rm_get_kernel # create reduced lda model from simulation data
I.silence_stdout # context manager and decorator to silence functions
I.cache # decorator to cache functions
I.cluster() # get distributed.Client() object for parallel execution
I.np # numpy
I.pd # pandas
I.dask # dask
I.distributed # distributed
I.sns # seaborn
# ...
Use autocompletion of IPython to find the other modules. To view the documentation, use a questionamrk, e.g.
I.mdb_init_simrun_general.init?
Quickstart: How to run a simulation?¶
To run single cell simulations, we first need a hoc morphology file of that cell. Additionally, we need to know, how to distribute ion channels on that morphology, such that the cell has the desired biophysical properties.
In the getting started subfolder of this repository, there is such a hoc-morphology and the corresponding parameter file.
- anatomical_constraints/86_L5_CDK20041214_nr3L5B_dend_PC_neuron_transform_registered_C2center.hoc
- biophysical_constraints/86_CDK_20041214_BAC_run5_soma_Hay2013_C2center_apic_rec.param
Let's have a closer look at the parameters.
from getting_started import getting_started_dir # path to getting started folder
cell_param = I.os.path.join(getting_started_dir, \
'biophysical_constraints', \
'86_CDK_20041214_BAC_run5_soma_Hay2013_C2center_apic_rec.param')
cell_param = I.scp.build_parameters(cell_param) # this is the main method to load in parameterfiles
cell_param is a nested dictionary. The biophysical parameters are described in the neuron key. Here, we specify ion channel conductances for all structures of the cell:
cell_param.neuron.keys()
The key filename points to the hoc morphology file. The others define ion channel conductances for the respective part of the cell. Let's figure out, what channels can be found in the basal dendrite:
cell_param.neuron.Dendrite
Here you see that we only have the passive conductance and If channels, both uniformly distributed along the dendrite. There are much more channels in the apical dendrite. A comprehensive description of the other parts of the cell parameter file can be found here: #todo
Now, we can do some current injection experiments to analyze the biophysical properties of the model.¶
# set up cell. This creates a lot of (interesting) diagnostic output. Remove the context manager to see it
with I.silence_stdout:
cell = I.scp.create_cell(cell_param.neuron)
Let's figure out some general properties of the cell
print 'total length = {:.0f} micron'.format(sum([sec.L for sec in cell.sections]))
print 'total dendritic length = {:.0f} micron'.format(sum([sec.L for sec in cell.sections \
if sec.label in ['Soma', 'Dendrite', 'ApicalDendrite']]))
print 'soma area = {:.0f} micron^2'.format(sum([sec.area for sec in cell.sections if sec.label == 'Soma']))
print 'apical dendrite area = {:.0f} micron^2'.format(sum([sec.area for sec in cell.sections \
if sec.label == 'ApicalDendrite']))
# load scaled hoc morphology
cell_param.neuron.filename = I.os.path.join(getting_started_dir, 'anatomical_constraints', \
'86_L5_CDK20041214_nr3L5B_dend_PC_neuron_transform_registered_C2center_scaled_diameters.hoc')
with I.silence_stdout:
cell = I.scp.create_cell(cell_param.neuron)
print 'apical dendrite area = {:.0f} micron^2'.format(sum([sec.area for sec in cell.sections \
if sec.label == 'ApicalDendrite']))
Now, we put a pipette at the soma of the cell and inject a short rectangular current
import neuron
h = neuron.h
iclamp = h.IClamp(0.5, sec=cell.soma)
iclamp.delay = 150 # give the cell time to reach steady state
iclamp.dur = 5 # 5ms rectangular pulse
iclamp.amp = 1.9 # 1.9 ?? todo ampere
%time I.scp.init_neuron_run(cell_param.sim, vardt=True) # run the simulation
%matplotlib inline
I.sns.set_style('ticks')
I.plt.plot(cell.tVec, cell.soma.recVList[0], c = 'k')
I.sns.despine()
How does the cell respond to different amplitudes of the step current?
cell_param.sim.tStop = 3000
iclamp.dur = 2000
iclamp.delay = 500
simresult = {}
for amp in [0.619, 0.793, 1.507]:
iclamp.amp = amp
I.scp.init_neuron_run(cell_param.sim, vardt=True)
simresult[amp] = I.np.array(cell.tVec), I.np.array(cell.soma.recVList[0])
fig = I.plt.figure(figsize = (9,3), dpi = 200)
for lv, k in enumerate(sorted(simresult.keys())):
ax = fig.add_subplot(1,3,lv+1)
ax.plot(*simresult[k], label = k, c = 'k')
ax.legend(loc = 'upper right')
ax.set_ylim([-90, 40])
ax.set_xlabel('t / ms')
if lv == 0:
ax.set_ylabel('Vm / mV')
I.sns.despine()
Now let's place a pipette on the apical trunk at a soma distance of 620 microns and inject current synchronously at this position and at the soma. For the current at the soma, we choose a rectangular pulse again. The current at the apical injection site should resemble the shape of a evoked post synaptic potential (epsp). Have a look at the folder mechanisms/synapses to see how epsp is implemented in hoc code.
# map between sections and somadistance
sec_dist_dict = {cell.distance_to_soma(sec, 1.0): sec for sec in cell.sections}
def get_section_at_distance(cell, dist):
dummy = {k - dist: v for k,v in sec_dist_dict.iteritems() if k > dist}
closest_sec = dummy[min(dummy)]
x = (dist - cell.distance_to_soma(closest_sec, 0.0)) / closest_sec.L
return x, closest_sec
cell_param.sim.tStop = 500
x, sec = get_section_at_distance(cell, 620)
iclamp = h.IClamp(0.5, sec=cell.soma)
iclamp.delay = 295 # give the cell time to reach steady state
iclamp.dur = 5 # 5ms rectangular pulse
iclamp.amp = 1.9 # 1.9 ?? todo ampere
iclamp2 = h.epsp(x, sec=sec)
iclamp2.onset = 300
iclamp2.imax = 0.5
iclamp2.tau0 = 1.0 # rise time constant
iclamp2.tau1 = 5 # decay time constant
cell.iclamp2 = iclamp2
%time I.scp.init_neuron_run(cell_param.sim, vardt=True) # run the simulation
I.plt.plot(cell.tVec, cell.soma.recVList[0], c = 'k')
I.sns.despine()
So far, we have only visualized the somatic membrane potential. Now, we create a short animation, where the whole dendrite is visualized (x-axis: soma distance of dendrite segment, y-axis: Vm)
mdb = I.ModelDataBase('/home/abast/my_mdb/') # create a ModelDataBase to save .png-frames of the animation
see section ModelDataBase for a more detailed descritpion
from IPython import display
if not 'burst_trail_video' in mdb.keys():
I.scp.init_neuron_run(cell_param.sim, vardt=False) # fixed step size necessary for animation
I.cell_to_animation(cell, \
xlim = [0,1500], \
ylim = [-90, 50], \
tstart = 295, \
tend = 350, \
tstep = 0.2, \
outdir = mdb.create_managed_folder('burst_trail_video'))
I.display_animation(mdb['burst_trail_video'].join('*', '*.png'), embedded=True)
We can also visualize current through ion channels. Let's have a look at the calcium current through low voltage activated calcium channels:
if not 'burst_trail_ca_current_video' in mdb.keys():
cell.record_range_var('Ca_LVAst.ica') # have a look at the mechanisms folder to find available range vars
I.scp.init_neuron_run(cell_param.sim, vardt=False)
I.cell_to_animation(cell, \
xlim = [0,1500], \
ylim = [-0.3, 0.1], \
tstart = 295, \
tend = 350, \
tstep = 0.2, \
outdir = mdb.create_managed_folder('burst_trail_ca_current_video'),\
range_vars = 'Ca_LVAst.ica')
I.display_animation(mdb['burst_trail_ca_current_video'].join('*', '*png'), embedded = True)
See section display_animations for caveats of creating animations
Step by step: Simulations with well-constrained synaptic input¶
Here, we place a cell in barrel cortex and compute an anatomical realization, which specifies, where synapses are located. Afterwards, we activate the synapses based on experimental measurements. Next, we import the simulation result in a ModelDataBase to analyze it. Whenever necessary, dask and distributed is used to parallelize the computation.
Step 1: Register cell morphology in the BC reference frame¶
From external resources (e.g. the NeuroMorph pipeline), you need a hoc-file with the dendrite morphology you want to place in the barrel cortex model. The coordinates in the hoc morphology file needs to be choosen according to the coordinate system of the barrel cortex model, i.e. the cell should sit at the desired location.
There is a morphology file in the getting_started/anatomical_constraints subfolder, which is registered such that the cell is in the center of the C2 column. For the following steps, this morphology is used:
I.os.listdir(I.os.path.join(getting_started_dir, 'anatomical_constraints'))
For building the anatomical model, we use '86_L5_CDK20041214_nr3L5B_dend_PC_neuron_transform_registered_C2center.hoc', for the simulation of the evoked activity, we need to use the morphology file with scaled apical trunk. Why the scaling is necessary is described above.
path_to_hoc = I.os.path.join(getting_started_dir, 'anatomical_constraints', \
'86_L5_CDK20041214_nr3L5B_dend_PC_neuron_transform_registered_C2center.hoc')
path_to_scaled_hoc = I.os.path.join(getting_started_dir, 'anatomical_constraints', \
'86_L5_CDK20041214_nr3L5B_dend_PC_neuron_transform_registered_C2center_scaled_diameters.hoc')
We copy both files to our ModelDataBase.
mdb = I.ModelDataBase('/home/abast/my_mdb/')
if not 'anatomical_constraints' in mdb.keys():
mdb.create_managed_folder('anatomical_constraints')
I.shutil.copy(path_to_hoc, mdb['anatomical_constraints'])
I.shutil.copy(path_to_scaled_hoc, mdb['anatomical_constraints'])
Sidenote: The copying step is optional, if we just want to run the simulation. To make the database self-containing (i.e. after copying it to another machine to an arbitrary path, all paths in all parameter files are stil valid and simulations can be rerun immediately), it is required that all files are in a ModelDataBase: In this case, we could use the following path to refer to the hoc file:
mdb://2017-10-19_3728_uOGyyhY/anatomical_constraints/some_file.hoc*
mdb.get_id()
mdb paths are described in more detail in section todo
Step 2: calculate positions of anatomical synapses with respect to the registered morphology¶
We use the singlecell_input_mapper to create an anatomical model of how that cell is integrated in barrel cortex. The following code calculates the position of synapses onto the example cell. Under the hood, 50 anatomical realizations will be computed. Out of these, the one that is closest to the average is choosen, which then can be refered to as "representative realization". The result is saved in the same folder as the hoc morphology. This takes about 4 hrs to compute, but you can continue with a precomputed result.
celltype = 'L5tt'
path = mdb['anatomical_constraints'].join('86_L5_CDK20041214_nr3L5B_dend_PC_neuron_transform_registered_C2center.hoc')
with I.silence_stdout: # this creates > 1 GB of diagnostic information
I.map_singlecell_inputs(path, celltype) # this takes some time
To use the precomputed result, just copy it from the getting started folder to the ModelDataBase:
from distutils.dir_util import copy_tree
path_to_anatomical_model = I.os.path.join(getting_started_dir, 'anatomical_constraints', \
'86_L5_CDK20041214_nr3L5B_dend_PC_neuron_transform_registered_C2center_synapses_20150504-1611_10389/')
silent = copy_tree(path_to_anatomical_model, mdb['anatomical_constraints'])
# adapt the path, if you have generated a new anatomical model
path_to_anatomical_model = mdb['anatomical_constraints'].join('86_L5_CDK20041214_nr3L5B_dend_PC_neuron_transform_registered_C2center_synapses_20150504-1611_10389')
What does the result of the map_singlecell_inputs script look like?
In the directory generated by the singlecell_input_mapper, there are the following files:
I.os.listdir(path_to_anatomical_model)
The most important files are the .con file and the .syn. These files are the relevant output of the SingleCellMapper for simulations of evoked activity.
The summary.csv file contains summary statistics about the connectivity between all presynaptic populations and our neuron.
The con file maps between presynaptic cells and synapses:
con_file_path = I.glob.glob(path_to_anatomical_model.join('*.con'))[0]
with open(con_file_path) as f:
print f.read()[:300]
The .syn file specifies the exact position of each synapse on the hoc morphology:
syn_file_path = I.glob.glob(path_to_anatomical_model.join('*.syn'))[0]
with open(syn_file_path) as f:
print f.read()[:300]
Here, section referes to the ID od the section in the cell object. x specifies, where along that section the synapse is placed. If x is 0, this is the beginning of the section, if x is one, this is the end of the section.
The subfolders in the directory contain Amira landmark files:
- for the location of presynaptic_somata by celltype in the folder
presynaptic_somata - for the location of synapses in the folders
total_synapses,soma_synapses,basal_synapses,apical_synapses
How to parallelize the generation of anatomical models?
We can use a distributed Cluster:
I.cluster()
# make the single cell mapper a delayed function
delayed_map_singlecell_inputs = I.dask.delayed(I.map_singlecell_inputs)
# call it with the morphologies
delayeds = [delayed_map_singlecell_inputs(p, 'L5tt') for p in morphology_paths]
# bundle everything in one delayed object
delayeds = dask.delayed(delayeds) # just do it ...
# compute the result, in this case, all cores on the local machine are used.
futures = I.cluster().compute(delayeds)
# visualize the progress
I.distributed.progress(futures)
Sidenotes
I.cluster() returns a distributed.Client object, which can be used to execute a dask graph using all cores on the local machine. I.cluster() also "prepares" the cluster such that it fits the needs of doing single cell simulations:
- Work stealing is disabled, because otherwise, tasks might run more than once. This should change in the near future, see this github issue. As soon as work-stealing is transactional, just comment out the respective lines in Interface.cluster, which disable work stealing.
- It imports matplotlib on all workers and sets the backend to 'agg'. This can be neccessary on compute servers that do nt provide the necessary graphical libraries, which would crash otherwise. After this, you can create figures on the computeservers.
- If I.cluster() is called more than once, it does not generate a new distributed.Client() object every time. Instead, it creates it when it is first called and returns the same client on all consecutive calls.
- You can set up a "real" cluster spanning hundrets of compute servers as described here. In this case, you would connect to the cluster in the following way:
I.cluster('server_name_on_which_scheduler_is_running:8786')
Step 3: define synapse activation¶
Now we need to activate the synapses according to experimental data, such that the synapse activation represents a passive touch experiment.
To do this, we need the following information
- A parameter file specifying charakteristics of the synaspses by celltype
- A parameterfile specifying the ongoing activity by celltype
- Parameterfiles specifying evoked response by stimulus, celltype and celllocation
Synapse characteristics and ongoing activity¶
Synapse characteristics and ongoing activity are defined here:
ongoing_template_param_name = I.os.path.join(getting_started_dir, 'functional_constraints/ongoing_activity/ongoing_activity_celltype_template_exc_conductances_fitted.param')
In this parameter file, the following keys are defined
ongoing_template_param = I.scp.build_parameters(ongoing_template_param_name)
ongoing_template_param.keys()
The most relevant information is specified in the network key:
print ongoing_template_param.network.keys()
Here, parameters are defined for each presynaptic celltype:
ongoing_template_param.network.L5tt
- interval: mean ongoing interspike interval
- releaseProb: chance that synapse gets activated if presynaptic cell gets activated
- synapses --> receptors --> glutamate_syn --> parameter:
- parameters for the NEURON mechanism defined here: mechanisms/channels/netglutamate.mod
- tau1: NMDA decay time
- tau2: NMDA rise time
- tau3: AMPA decay time
- tau4: AMPA rise time
- parameters for the NEURON mechanism defined here: mechanisms/channels/netglutamate.mod
- synapses --> receptors --> weight
- maximum conductance of the synapse for AMPA and NMDA, respectively
PSTHs of evoked activity for passive touch by celltype and celllocation¶
Parameterfiles can be found here:
evokedPrefix = I.os.path.join(getting_started_dir, 'functional_constraints/evoked_activity/PW_SuW_RF_CDK/')
excitatory_PSTHs = [fname for fname in I.os.listdir(evokedPrefix) if fname.endswith('PSTH_UpState.param')]
inhibitory_PSTHs = [fname for fname in I.os.listdir(evokedPrefix) if fname.endswith('active_timing_normalized_PW_1.0_SuW_0.5.param')]
excitatory_PSTHs
Let's have a look at a specific parameterfile:
example_PSTH = I.scp.build_parameters(I.os.path.join(evokedPrefix, excitatory_PSTHs[-5]))
example_PSTH.keys()
We have entries for each homecolumn, each containing the respective evoked PSTH in a C2 stimulus scenario:
example_PSTH.L5tt_C2.keys()
for k,v in example_PSTH.L5tt_C2.iteritems():
print '{}: {}\n'.format(k,v)
Creating a parameterfile that fully defines a specific stimulus¶
Next, we combine all data, to one parameterfile, which then contains all information to describe the activatio of presynaptic cells and synapses during a passive whisker touch scenario, given our cellmorphology and the anatomical network realization. To so, we use the create_evoked_network_parameter function:
# the whisker stimulus we want to simulate.
whisker = 'C2' # Cell is in C2 --> this is a principal whisker stimulus scenario
# defined above
ongoing_template_param_name
syn_file_path
con_file_path
# cell number spreadsheet generated in anatomical realization step.
cell_number_file_name = path_to_anatomical_model.join('NumberOfConnectedCells.csv')
# output_path
if not 'network_param' in mdb.keys(): mdb.create_managed_folder('network_param')
out_file_name = mdb['network_param'].join('C2_stim.param')
I.create_evoked_network_parameter(ongoing_template_param_name, cell_number_file_name, \
syn_file_path, con_file_path, whisker, out_file_name)
Let's also generate parameterfiles for surround whisker stimuli!
with I.silence_stdout:
for whisker in ['B1', 'B2', 'B3', 'C1', 'C3', 'D1', 'D2', 'D3']:
I.create_evoked_network_parameter(ongoing_template_param_name, cell_number_file_name, \
syn_file_path, con_file_path, whisker,\
mdb['network_param'].join('{}_stim.param'.format(whisker)))
what does the result look like?¶
network_param = I.scp.build_parameters(mdb['network_param'].join('C2_stim.param'))
network_param.keys()
print network_param.network.keys()[:10]
We just printed a subset of the keys here. Lets have a look at what the parameterfile specifies about L5tt cells in C2
for k, v in network_param.network.L5tt_C2.iteritems():
print '{}: {}\n'.format(k,v)
As you see, all information is integrated
- network --> [population] --> celltype --> pointcell --> offset: timepoint at which the evoked activity starts
Step 4: Determine channel densities such that the cell exhibits certain biophysical properties¶
Currently, the parameters are determined using a MOO approach. The result is writen manually in the cell parameter file, which we have already seen above.
# We copy it to the model data base
if not 'biophysical_constraints' in mdb.keys():
cell_parameter_file = I.os.path.join(getting_started_dir, 'biophysical_constraints', \
'86_CDK_20041214_BAC_run5_soma_Hay2013_C2center_apic_rec.param')
I.shutil.copy(cell_parameter_file, mdb.create_managed_folder('biophysical_constraints'))
cell_parameter_file = mdb['biophysical_constraints'].join('86_CDK_20041214_BAC_run5_soma_Hay2013_C2center_apic_rec.param')
Step 5: run the simulation¶
we simulate one single trail for each whisker stimulus and save the result in the database under the simrun key:
delayeds = []
if 'simrun' in mdb.keys():
del mdb['simrun']
mdb.create_managed_folder('simrun')
for whisker in ['B1', 'B2', 'B3', 'C1', 'C2', 'C3', 'D1', 'D2', 'D3']:
network_file = mdb['network_param'].join('{}_stim.param'.format(whisker))
cell_parameter_file # defined above
dir_prefix = mdb['simrun'].join('stim_{}'.format(whisker))
nSweeps = 1 # number of consecutive simulations per process
nprocs = 1 # number of processes simulating that task in parallel
tStop = 300 # stop the simulation at 300ms
d = I.simrun_run_new_simulations(cell_parameter_file, network_file, dirPrefix = dir_prefix, \
nSweeps = nSweeps, nprocs = nprocs, tStop = tStop, silent = False)
delayeds.append(d)
I.simrun_run_new_simulations returns dask.delayed objects, i.e. they result describes how the computation can be run, but it is not computed yet. We send it to our local cluster:
futures = I.cluster().compute(delayeds)
The computation happens in the background and does not block IPython. We can visualize the progress:
I.distributed.progress(futures)
Now, the simrun folder in the model database contains the simulation results
I.os.listdir(mdb['simrun'])
Step 6: load computation results to a ModelDataBase for further analysis¶
This gives convenient access to the data using the widely used pandas and dask libraries. It also converts the data to an optimized binary format (msgpack with blosc compression, also the synapse and cell activation data is categorized, i.e. repeated values are replaced by an integer). This gives the following advantages:
- reduced space on harddisk (25% of original data)
- fast access (... s to load in the synapse activation of 1000 simulation trails compared to .. It also provides
I.mdb_init_simrun_general.init(mdb, mdb['simrun'], get = I.cluster().get)
I.mdb_init_simrun_general.optimize(mdb, get = I.cluster().get)
mdb.keys()
Now, the following keys are set:
- mdb['synapse_activation'] activation of synapses for each simulationtrail, type: dask:DataFrame
- mdb['cell_activation'], activation of presynaptic cells for each simulationtrail, type: dask:DataFrame
- mdb['voltage_traces'], soma voltage traces, type: dask:DataFrame
- mdb['spike_times'], spike times, type: pandas.DataFrame
- mdb['dendritic_recordings'], voltage traces on the dendrite as specified in the cell parameter file, type: ModelDataBase containing a seperate dask.DataFrame for every recording location
- mdb['metadata'], original filenames, type: pandas.DataFrame
- mdb['parameterfiles'], mapping between simulation trial and used parameterfiles, type: pandas.DataFrame
- mdb['parameterfiles_cell_folder'], folder containing the cell parametefiles, type: ManagedFolder
- mdb['parameterfiles_network_folder'], folder containing the network parameterfiles, type: ManagedFolder
- mdb['sim_trail_index'], all simulation trials, type: pandas.Series
- mdb['simresult_path'], original path from which the simulation has been loaded, type: str
Let's have a look at the synapse activation data frame:
mdb['synapse_activation']
This is a dask dataframe object, which points to the actual data which is stored in the database using the msgpack binary format. The data is split up into chunks of reasonable file size. At this timepoint, no data has been loaded. We can view the first 2 of the dataset, which only loads the first chunk of data:
mdb['synapse_activation'].head(2)
We can select data using the index:
mdb['synapse_activation'].loc['stim_B1/results/20171113-1614_seed1467457903/000000']
This returns a dask dataframe again, i.e. the query is not computed yet. To do this, call compute:
# this prints out a table containing all synapse activations during the specified trail
mdb['synapse_activation'].loc['stim_B1/results/20171113-1614_seed1467457903/000000'].compute()
What is the average distance of synapses from the soma by celltype?
# load dask dataframe
sa = mdb['synapse_activation']
# assign a new column: the celltype is the first part of the synapse_type column
sa['celltype'] = sa.map_partitions(lambda x: x.synapse_type.str.split('_').str[0])
# groupb by this value, compute the mean soma distance in each group, compute the result on the cluster
sa.groupby('celltype').apply(lambda x: x.soma_distance.mean()).compute(get = I.cluster().get)
It is easy to create an animation of a simulation trail from the database:
if not 'D3_stim_animation' in mdb.keys():
cell = I.simrun_simtrail_to_cell_object(mdb, 'stim_D3/results/20171113-1614_seed919944399/000000')
I.cell_to_animation(cell, \
xlim = [0,1500], \
ylim = [-90, 50], \
tstart = 245-25, \
tend = 245+25, \
tstep = 0.2, \
outdir = mdb.create_managed_folder('D3_stim_animation'))
I.display_animation(mdb['D3_stim_animation'].join('*', '*png'), embedded = True)